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Abstract  The  present  study  is  concerned  with  the  numerical  simulation  of  high  Reynolds 
number  weakly  separated  steady  laminar  flows.  It  is  shown  that  a  "well  designed''  code  can  solve 
the  complete  Navier-Stokes  (NS)  equations  and  the  "more  suitable?"  parabolized  Navier-Stokes 
(PNS)  equations  with  the  same  convergence  rate,  so  that  solving  the  full  NS  equations  is 
recommended  when  dealing  with  a  new  problem.  Furthermore,  solutions  to  the  classical 
boundary  layer  equations  in  their  vorticity-stream  function  form  are  obtained,  which  are  regular 
thru  the  separation  point,  i.  e.,  do  not  encounter  the  Goldstein  singularity  at  separation.  Finally, 
for  a  very  typical  high  Reynolds  number  weakly  separated  flow,  it  is  shown  that,  in  the  presence 
of  non  negligible  skewness  in  the  body  oriented  computational  grid,  a  PNS-type  approximation 
still  provides  very  reliable  solutions,  whereas  an  interacting  boundary  layer-type  model  is  plagued 
by  severe  errors. 
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Introduction 

High  Re  (Reynolds  number)  attached  flows  have  been  computed  routinely  for  several  decades 
using  the  classical  boundary  layer  (BL)  theory  of  Prandtl:  the  inviscid  flow  past  the  body  of 
interest  is  computed  at  first  and  then  the  BL  equations  are  solved,  using  the  inviscid  pressure 
gradient  as  a  datum,  to  evaluate  the  skin  friction  and  thus  the  aerodynamic  drag.  Eventually, 
the  displacement  thickness  obtained  by  the  BL  equations  is  added  to  the  body  to  provide  an 
"augmented  body"  around  which  the  inviscid  flow  is  recomputed  and  so  on,  according  to  a 
rigorous  asymptotic  expansion  theory  provided  by  Van  Dyke  [l]. 

Unfortunately,  the  classical  BL  theory  breaks  down  in  the  presence  of  flow  separation,  so  that 
a  lot  of  computational  effort  has  been  devoted  to  the  development  of  numerical  methods  for 
solving  the  Navier-Stokes  (NS)  equations.  However,  theoretical  studies  (see,  eg,  [2])  as  well  as 
numerical  investigations  (see,  eg,  [3,  4])  have  shown  that,  for  high  Re  weakly  separated  flows,  the 
terms  contained  in  the  classical  BL  equations  are  sufficient  to  model  the  flow  field  even  in  the 
separated  region.  It  is  now  understood  that  the  failure  of  the  classical  BL  theory  (when  the 
pressure  gradient  is  prescribed)  is  due  to  the  lack  of  upstream  propagation  required  to  account  for 
the  "elliptic”  nature  of  the  separation  phenomenon;  and  in  fact,  solutions  to  the  BL  equations, 
which  are  regular  thru  the  separation  point,  have  been  obtained  by  prescribing  either  the 
displacement  thickness  [5],  or  the  wall  shear  [6],  in  place  of  the  pressure  gradient. 

More  in  general,  useful  solutions  for  high  Re  separated  flows  can  be  obtained  from  an 
appropriate  subset  of  the  NS  equations,  provided  that  the  "elliptic”  nature  of  the  strong  viscous 
inviscid  interaction  phenomenon  is  correctly  accounted  for.  Among  the  many  such  approximate 
equations  used  in  the  literature,  the  most  widely  employed  and  investigated  appear  to  be  the  so- 
called  parabolised  Navier-Stokes  (PNS)  and  interacting  boundaiy  layer  (IBL)  equations.  Starting 
from  the  NS  equations  written  in  a  body  oriented  orthogonal  coordinate  system  the  PNS 
equations  are  easily  obtained  by  dropping  the  diffusion  terms  in  the  streamwise  direction.  The 
IBL  equations  are  obtained  instead  by  coupling  the  classical  BL  equations  with  a  pressure 
interaction  law,  which  accounts  for  the  variations  in  the  inviscid  (outer)  pressure  field  induced  by 
the  viscous  effects.  The  PNS  equations  can  then  be  considered  from  an  asymptotic  point  of  view 
as  a  composite  set  of  equations  uniformly  valid  in  both  the  inner  (viscous)  and  the  outer  (inviscid) 
regions  and  therefore  "equivalent”  to  the  complete  NS  equations.  The  IBL  equations  are  valid 
instead  only  in  the  viscous  region,  whereas  the  inviscidTlow  is  computed  by  means  of  another  set 
of  equations;  and  the  viscous-inviscid  interaction  phenomenon  is  accounted  for  iteratively,  e.  g., 
by  solving  the  inviscid  and  viscous  flow  equations  using  the  same  displacement  thickness 
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(prescribed  appropriately)  and  repeating  the  calculations  until  both  sets  of  equations  provide  the 
same  pressure  distribution  along  the  displacement  body.  For  the  case  of  the  vorticity-stream 
function  equations,  of  interest  here,  the  PNS  equations  are  obtained  by  simply  dropping  the 
stream  wise  diffusion  term  in  the  vorticity  equation,  so  that  they  are  also  known  in  the  literature 
as  the  parabolized  vorticity  (PV)  equations  [3,  4].  The  classical  BL  equations  are  obtained  by 
dropping  the  streamwise  derivative(s)  also  in  the  stream  function  equation,  according  to  the 
asymptotic  order  of  magnitude  analysis  of  Prandti,  and  the  IBL  equations  are  therefore  obtained 
by  coupling  the  BL  equations  in  the  viscous  region  (where  the  vorticity  is  nonzero)  with  the 
Laplace  equation  for  the  stream  function  in  the  inviscid  region,  so  as  to  allow  for  the  viscous- 
inviscid  interaction  phenomenon. 

To  date,  a  large  number  of  numerical  studies  have  shown  that  for  many  flow  fields  of  practical 
interest,  namely  for  all  of  the  high  Re  weakly  separated  flows,  the  PNS  or  even  the  IBL  equations 
provide  a  satisfactory  answer,  so  that  it  would  appear  as  an  unnecessary  waste  of  computational 
effort  to  resort  to  the  complete  NS  equations,  see,  e.  g.,  [7-10).  Furthermore,  due  to  the 
pioneering  theoretical  work  of  Smith  and  Chen  [11-13],  it  starts  to  appear  more  and  more  evident 
that  the  IBL  equations  are  capable  of  predicting  also  the  massive  separation  phenomenon  and  the 
stall,  as  it  has  also  been  shown  numerically  by  Rothmayer  and  Davis  [14]. 

However,  the  NS  equations  still  remain  the  only  numerical  tool  capable  of  verifying  the 
correctness  of  the  new  asymptotic  theories  and,  if  solved  properly,  they  require  a  computational 
effort  which  is  of  the  same  order  of  magnitude  as  that  required  by  the  PNS,  if  not  by  the  IBL 
equations.  In  fact,  although  obvious,  it  does  not  seem  to  be  generally  understood  that  a 
numerical  method  for  solving  the  NS  equations  for  high  Re  weakly  separated  flows  has  to  be  very 
similar  in  nature  to  efficient  PNS  or  IBL  codes;  i.  e.,  it  has  to  account  for  downstream- 
propagation  and  transversal-viscous  effects  implicitly  and  only  for  the  upstream-influence  in  an 
iterative  fashion.  For  such  a  code,  the  convergence  rate  will  be  the  same  (see,  e.  g.,  [10])  whether 
only  the  terms  present  in  the  PNS  (or  even  the  IBL)  equations  are  computed  or  all  of  the  terms 
in  the  complete  NS  equations  are  accounted  for;  and  the  extra  computational  cost  required  to 
compute  the  negligible  streamwise  diffusion  terms  will  be  insignificant  with  respect  to  that 
required  to  solve  the  PNS  equations.  But  what,  if  one  of  these  term  is  essential  to  model  a 
locally  important  phenomenon?  It  is  worth  pointing  out  that  many  IBL  codes  only  solve  the  BL 
equations  (iteratively)  in  the  viscous  region,  with  a  far  Held  condition,  obtained  by  computing  the 
inviscid  flow  past  the  displacement  body  according  to  thin  airfoil  theory,  accounting  for  the 
interaction  phenomenon  (see,  e.  g.,  [14]).  These  methods,  of  course,  may  require  significantly  less 
computational  effort  than  codes  solving  the  complete  NS  equations,  especially  for  the  case  of 
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supersonic  outer  flow. 

The  aim  of  this  work  is  to  demonstrate  that  for  high  Re  weakly  separated  flows,  a  well 
designed  NS-PNS  code  can  compute  the  NS  or  the  PNS  equations  with  exactly  the  same 
convergence  rate  and  at  about  the  same  computational  cost,  so  that,  when  solving  a  novel 
problem,  the  additional  cost  of  carrying  all  of  the  terms  contained  in  the  NS  equations  is  minimal 
and  certainly  justified  by  the  fact  that  the  model  employed  accounts  for  all  meaningful 
phenomena.  Of  course,  by  no  means  it  is  implied  that  solving  the  NS  equations  by  "brute  force" 
is  an  appropriate  investigating  tool!  Asymptotic  analyses  are  essentia]  indeed,  to  prescribe 
appropriate  computational  grids  and  the  location  and  type  of  the  farfield  boundary  conditions,  as 
well  as  to  validate  (and  be  validated  by)  numerical  solutions. 

Furthermore,  the  present  study  aims  to  demonstrate  that,  when  solving  the  vorticity  stream 
function  equations,  solutions  for  high  Re  separated  flows,  which  are  regular  thru  the  separation 
point,  can  be  obtained  by  solving,  in  the  viscous  region,  simply  the  classical  BL  equations, 
namely,  by  dropping  the  streamwise  derivatives  also  in  the  stream  function  equation.  In 
particular,  for  the  case  of  the  high  Re  flow  in  a  symmetric  channel  with  a  smooth  expansion 
proposed  by  Roache  [15],  which  is  of  parallel  type,  namely  has  negligible  streamwise  curvature 
effects,  these  solutions  of  the  BL  equations  will  be  shown  to  be  in  very  good  agreement  with 
results  of  the  PNS  and  of  the  complete  NS  equations. 

Finally,  for  a  very  typical  high  Re  weakly  separated  flow,  already  studied  by  several  authors 
[7-10],  in  the  presence  of  skewness  of  the  body  fitted  coordinate  system,  it  is  shown  that  PNS- 
type  equations  still  provide  results  in  good  agreement  with  the  complete  NS  solutions,  whereas 
IBL-type  equations  experience  marked  errors  due  to  the  incorrect  representation  of  the 
streamwise  pressure  gradient. 

Governing  Equations  and  Numerical  Technique 

The  vorticity-stream  function  equations  are  given  in  a  system  of  body-oriented  curvilinear 
coordinates  as: 

'  2^(f,  +  +  <%  +  Wg)/R»  =  0  (1) 

QV« '  +  7i’™  +  ^  =o 


(2) 
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In  eqns  (1-2)  -u>  i«  the  vorticitv,  t!'  is  the  stream  function,  £  and  17  are  the  streamwise  and 
transversal  coordinates,  t  is  the  time,  subscripts  indicate  partial  derivatives  and  J,  a,  0,  7,  a  and 
r  are  the  jacobian  and  the  scale  factors  of  the  mapping  of  the  physical  coordinates  x,  y  into  the 
computational  coordinates  £,  t)  [16].  Furthermore,  the  terms  underlined  are  those  which  are  to  be 
dropped  to  provide  the  PNS  equations  and  those  underlined  with  a  broken  line  are  the  additional 
terms  which  also  need  to  be  dropped  to  obtain  the  classical  BL  equations.  Of  course,  BL  theory 
strictly  requires  the  use  of  orthogonal  coordinates,  for  which  0  in  eqns  (1-2)  is  identically  zero. 
However,  the  PNS  equations  have  been  shown  to  be  valid  also  in  a  "Mildly  Nonorthogonal 
Coordinate  System”  [10],  so  that  the  aforementioned  definition  of  PNS  and  BL  equations  is  used 
here.  Eqn  (2),  with  the  vorticity  term  eliminated,  obviously  is  the  Laplace  equation  for  the 

stream  function  governing  the  inviscid  irrotational  portion  of  the  flow  field. 

<\ 

In  the  present  study  the  incremental  block-line  Gauss-Seidel  method  proposed  in  Ref.  17  is 
used  as  an  efficient  numerical  tool  for  solving  eqns  (1-2)  numerically.  The  governing  equations  are 
discretized  and  linearized  in  time  using  a  two  level  implicit  Euler  scheme  and  the  delta  approach 
of  Beam  and  Warming  [18]  to  give: 

Au/Al  +  (V’^n  Alj^  +  -  u^Mt^)/ J 

-  +  (tAu ;  +  rAu>.)/Re  = 

*  K  V  w")*  -  U’^g/J  +  (au^n  •  20U;*  +  iwwn  +  tru*  +  Twfn)/Re  (3) 

Arl>/At  -  aAi'^  -  *  Au  = 

oV»ffn  -  2P4>(t)n  +  +  <^n  +  +  w®  (4) 

where  a  relaxation-like  time  derivative  has  been  added  to  the  stream  function  equation  to 
parabolize  it,  At  is  the  time  step  (which  can  be  different  in  eqns  (3)  and  (4)  and  can  also  vary  at 
every  iteration),  Au>  —  un+1  -  wn,  the  superscripts  n+1  and  n  indicating  the  new  and  old  time 
levels  tn+1  and  tn,  etc.  Notice  that  the  mixed  derivatives  are  evaluated  at  the  old  time  level  t“, 
i.  e.,  explicitly  and  therefore  do  not  appear  in  delta  form.  The  presence  of  such  terms  in  the 
vorticity  equation  of  the  complete  NS  equations  (when  a  nonorthogonal  computational  grid  is 
employed)  can  be  critical  insofar  as  they  may  reduce  the  stability  and  therefore  the  convergence 
rate  of  the  numerical  method. 

Eqns  (3-4)  are  then  discretized  in  space  using  second-order-accurate  central  differences 
throughout  -  except  for  the  incremental  convective  terms  of  the  vorticity  equation  in  the  left 


hand  side  of  eqn  (3),  which  are  approximated  using  first-order-accurate  upwind  differences  -  and 
solved  by  a  classical  line  Gauss-Seidel  relaxation  method,  marching  from  left  to  right  (see  [17],  for 
details).  It  is  noteworthy  that  by  using  a  deferred  corrector  strategy  for  the  advection  terms  in 
eqn  (3),  which  are  approximated  using  upwind  differences  in  the  LHS  and  central  differences  in 
the  RHS,  upstream  propagation  is  allowed,  although  explicitly,  even  when  solving  the  classical  BL 
equations,  i.  e.,  when  all  of  the  elliptic  terms  are  dropped  from  the  equations. 

Results 

The  first  problem  considered  in  this  study  is  the  weakly  separated  channel  flow  proposed  by 
Roache  [15].  The  geometry  of  the  channel,  together  with  the  appropriate  boundary  conditions 
and  the  orthogonal  41x21  grid  employed  in  this  study  is  given  in  Figure  1  for  the  case  Re  =*=  10. 
Roache  has  shown  that  if  the  length  of  the  channel  is  increased  proportionally  to  Re,  for  Re  >  > 
1,  the  solution  takes  on  a  quasi-self-similar  form;  i.  e.,  the  skin  friction  at  the  wall  becomes 
independent  of  Re,  when  plotted  versus  x/Re.  As  such,  the  present  problem  is  very  suitable  to 
assess  the  capability  of  a  numerical  technique  to  compute  high  Re  weakly  separated  flows;  more 
importantly,  for  Re  >>  1  ,  the  length  of  the  separation  bubble  is  found  to  increase 
proportionally  to  Re,  whereas  its  height  remains  constant,  so  that,  from  an  order  of  magnitude 
analysis,  it  can  be  argued  that: 


u  -  0(1) 

(5) 

v  *  0(l/Re) 

(6) 

d/dy  -  0(1) 

(7) 

d/dx  =  0(1  /Re) 

(8) 

where  0(  )  indicates  the  standard  asymptotic  order  of  magnitude  symbol.  If,  for  large  values 
of  Re,  equations  (5-8)  are  true,  than  the  PNS  equations  and  even  the  classical  BL  equations 
should  provide  results  practically  coincident  with  those  computed  using  the  complete  NS 
equations.  The  same  code  was  therefore  used  to  solve  the  three  aforementioned  sets  of  equations, 
for  the  case  Re  **  10®,  using  the  41x21  mesh  depicted  in  Figure  1.:  the  results  for  the  vorticity  at 
the  wall  were  found  to  be  identical  (to  machine  accuracy,  i.  e.,  6  significant  digits),  with  no 
indication  of  the  singularity  at  separation  appearing  when  solving  (iteratively)  the  classical  BL 
equations.  In  order  to  better  assess  this  last  very  important  issue,  as  well  as  the  accuracy  of  the 
present  numerical  results,  a  grid  refinement  study  was  conducted  in  both  the  x  (quasi-stream wise) 
and  the  y  (quasi- normal)  directions,  using  meshes  with  61x21,  81x21,  41x41  and  41x61  gridpoints, 


the  horizontal  mesh-size  being  always  uniform  and  the  normal  one  increasing  at  a  constant  rate 
equal  to  1.1,  starting  from  the  wall.  The  results  of  the  mesh  refinement  studies  are  given  in 
Figures  2  and  3  as  the  vorticity  distributions  along  the  wall:  it  clearly  appears  that  the  basic 
41x21-mesh  solution  is  already  "exact”,  to  plotting  accuracy.  Again,  the  NS,  PNS  and  BL  results 
coincided  to  machine  accuracy  in  all  cases  and  no  indication  of  singularity  at  separation  was  ever 
encountered.  It  is  noteworthy  that  for  the  present  flow  case,  the  length  of  the  channel  being 
equal  to  Re/3,  the  longitudinal  step  size  is  always  several  orders  of  magnitude  larger  than  the 
order-one  BL  thickness,  so  that  also  an  "ill- posed"  marching  method  could  possibly  recover  the 
correct  solution  [15],  by  swamping  out  all  branching  solutions  due  to  truncation  errors.  Therefore, 
the  Re  —  100  flow  case  was  also  considered  in  this  study:  the  BL  equations  were  solved  using  two 
meshes  having  81x21  and  161x21  gridpoints.  The  two  solutions  were  found  to  coincide  with  each 
other  and  (after  proper  scaling)  wi^>  those  of  Figures  2  and  3,  to  plotting  accuracy.  In 
conclusion,  the  BL  equations,  when  solved  in  an  iterative  way  which  allows  for  upstream 
propagation,  clearly  do  not  encounter  any  singularity  at  the  separation  point.  This  result  is  not 
surprising  insofar  as  the  value  of  the  stream  function  is  prescribed  as  a  boundary  condition  at  the 
centerline  of  the  channel;  such  a  condition  is  equivalent  tc  prescribing  the  displacement  thickness 
in  the  pressure-velocity  BL  equations,  which  has  been  shown  to  produce  separated  flow  solutions 
which  are  regular  thru  the  separation  point  [5],  Finally,  the  convergence  rate  for  all  three  sets  of 
equations  and  the  same  mesh  was  always  found  to  be  identical  and  is  provided  in  Figure  4,  where 
the  average  variation  of  the  vorticity  between  two  successive  iterations  is  plotted  versus  the 
iteration  number.  It  is  noteworthy  that  the  convergence  rate  of  the  method  is  remarkable  and 
just  about  independent  of  the  mesh  employed  for  the  calculations;  this  last  result  is  very  atypical 
of  any  relaxation  method,  so  that  it  is  most  likely  due  to  the  one-dimensional  (parallel)  nature  of 
the  flow. 

In  order  to  compare  the  present  NS  and  PNS  solvers  for  an  "extreme"  situation,  the  Re  ■=*  10 
flow  was  also  considered.  For  such  a  case,  the  geometry  is  considerably  distorted,  see  Figure  1, 
and  the  Reynolds  number  is  very  low,  so  that  it  is  reasonable  to  expect  large  differences  between 
the  NS  and  the  PNS  solutions,  whereas  the  BL  equations  are  obviously  meaningless.  However,  the 
present  NS  and  PNS  solutions  still  provide  wall  vorticity  distributions,  see  Figure  5,  which  are 
qualitatively  very  similar.  This  result  provides  further  evidence  of  the  broad  range  of 
applicability  of  the  PNS  equations. 

For  all  of  the  calculations  above,  an  appropriate  Couette  flow  was  used  at  every  longitudinal 
location  as  the  initial  condition  and  the  time  step  was  chosen  equal  to  Re,  initially,  and  then 
updated  after  every  iteration  by  dividing  its  initial  value  times  the  vorticity  error.  Therefore,  for 


the  particular  case  Re  =  106.  the  time  steps  were  so  large  that  the  time  derivatives  played  no 
role  in  the  computational  process. 

A  second  well  known  high  Re  weakly  separated  flow  problem  was  considered  in  the  present 
study,  namely  the  flow  in  a  symmetric  diffuser  already  studied  by  several  authors  [7-10],  in  order 
to  compare  the  NS  and  PNS  solutions  for  a  more  general  (non-parallel)  flow  and  in  the  presence 
of  local  skewness  of  the  computational  grid,  as  well  as  to  assess  the  efficiency  of  the  present  code 
with  respect  to  several  high  Re  weakly  separated  flow  solvers  [7-10].  Two  flow  cases  were 
considered,  whose  geometry  and  computational  grid  are  given  in  Figures  6a  and  6b.  In  both  cases 
the  nondimensional  semiheight  of  the  channel  is  equal  to  1,  the  Reynolds  number,  defined  with 
respect  to  such  a  reference  length  and  to  a  unitary  velocity,  is  equal  to  6250  and  the  inlet  flow  is 
a  Blasius  boundary  layer  produced  by  a  flat  plate  whose  leading  edge  is  two  units  before  the 
entrance  of  the  channel;  also,  the  midplane  of  the  channel  is  located  at  y  =  1,  whereas  the 
(bottom)  wall,  yw(x),  is  prescribed  as  follows: 

0  for  -1  <  x  <  0 

yw(x)  =  A[x“(3  -  2x )]  for  0  £  x  £  1  (9) 

A  for  1  <  x  <  xoul 

In  the  first  flow  case  A  =  -  0.08  and  xqu1  =  3,  whereas  in  the  second  more  severe  one,  A  — 
-  0.16  and  xQut  —  6.  The  computational  meshes  employed  in  this  study  contain  81  and  71 
equally  spaced  gridpoints  in  the  x  direction,  respectively,  and  41  gridpoints  in  the  y  direction, 
stretched  in  such  a  way  that  the  vertical  mesh  site  grows  at  a  constant  rate  of  1.08,  starting  from 
the  wall.  In  Figure  6a  only  every  other  gridline  is  plotted,  whereas  in  Figure  6b  all  of  the  vertical 
coordinate  lines  are  shown.  The  mesh  in  Figure  6a  was  chosen  fine  enough  to  provide  reasonably 
accurate  answers  to  compare  with  previously  published  results,  whereas  the  mesh  in  Figure  6b  is 
not  considered  adequate  to  provide  accurate  solutions,  but  valuable  to  provide  a  meaningful  one- 
to-one  comparison  between  NS  and  PNS  results  obtained  using  the  same  numerical  method  and 
mesh.  It  is  noteworthy  that  in  the  range  0  ^  x  <,  1,  the  mesh  is  nonorthogonal,  so  that  the 
mixed  derivatives  in  eqns  (1-4)  may  be  expected  to  play  a  considerable  role,  especially  for  the 
second  more  severe  flow  situation.  For  both  flow  cases,  results  were  obtained  using  the  same 
computer  code  to  solve  the  NS  and  the  PNS  equations.  The  initial  condition  was  taken  to  be  the 
inlet  boundary  condition  in  the  entire  flow  field.  Symmetry  was  enforced  by  prescribing  the 
value  of  the  stream  function  at  the  center  of  the  diffuser  (y  =*  1),  where  a  tero  vorticity  was  also 
prescribed.  At  the  wall  the  standard  no-slip  zero-injection  conditions  were  imposed  and,  finally, 
the  BL  equations  were  solved  at  the  outlet  station,  with  the  first  derivatives  evaluated  as  three- 
point  second-order-accurate  backward  differences  The  initial  time  step  was  taken  to  be  equal  to 
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0.1  for  the  vorticity  equation  and  1000  for  the  stream  function  equation.  After  every  iteration 
the  time  steps  were  adjusted  by  dividing  the  initial  values  by  the  average  vorticity  and  stream 
function  variations  (residuals),  respectively.  An  averrelaxation  factor  (8]  of  1.5  was  also  used,  for 
the  stream  function  equation  only.  All  of  these  values  are  the  first  and  oniy  ones  used,  so  that  no 
optimisation  of  the  convergence  rate  of  the  method  was  pursued.  For  the  First  flow  case,  IBL- 
type  equations  were  also  solved,  by  dropping  the  streamwise  derivatives  also  in  the  stream 
function  equation,  in  the  viscous  region,  while  retaining  it  in  the  outer  inviscid  flow  (namely,  in 
the  upper  half  of  the  computational  gridpoints). 

The  results  for  the  First  flow  case  are  given  in  Figures  7  and  8,  where  the  convergence  history 
for  the  NS  and  PNS  equations  and  the  vorticity  distributions  along  the  wall  for  all  of  the  three 
sets  of  equations  are  given,  respectively.  The  corresponding  results  for  the  second  flow  case  are 
given  in  Figures  9  and  10,  respectively.  From  the  results  in  Figures  7-10  the  following 
conclusions  emerge.  The  NS  and  PNS  equations  provide  practically  identical  results  and 
converge  at  the  same  rate.  The  IBL-type  equations  are  plagued,  instead,  by  local  errors  near  the 
curvature  discontinuities  at  the  wall,  where  the  pressure  gradient  induced  by  the  inviscid  turning 
of  the  streamlines  is  not  accounted  for  correctly.  It  is  noteworthy  that  for  the  second  flow  case 
the  IBL-type  equations  still  converge  without  difFiculties,  but  the  Final  solution  is  completely 
destroyed  by  spatial  oscillations  produced  by  the  more  severe  curvature  discontinuities.  Also,  the 
convergence  rate  of  the  IBL-type  equations,  not  reported  in  Figures  7  and  9,  is  slightly  slower 
than  that  of  the  NS  and  PNS  equations,  due  to  the  local  larger  gradients  to  be  computed. 

The  present  solution  given  in  Figure  8,  after  proper  scaling,  is  found  to  be  in  reasonable 
agreement  with  those  of  Edwards  and  Carter  [7)  and  of  Hoffman  [10],  but  predicts  a  slightly 
longer  separation  region.  The  PNS  equations  were  therefore  solved  using  a  Finer  121x01  mesh;  the 
results,  also  given  in  figure  8,  are  almost  identical  to  those  obtained  using  the  coarser  mesh.  The 
present  solutions  are  thus  considered  "correct",  to  plotting  accuracy,  and  the  slight  discrepancy 
with  the  aforementioned  solutions  is  believed  to  be  due  to  their  First-order- accuracy  in  the 
streamwise  direction.  For  the  second  flow  case,  these  authors  do  not  report  the  skin  friction 
distribution.  However  the  length  and  position  of  the  separated  region  can  be  compared  using  the 
streamlines  plot  of  Ref  7.  Again,  the  present  solutions  predict  a  separation  region  longer  than 
that  computed  by  Edwards  and  Carter  (7);  however,  in  this  case  the  present  mesh  is  not 
considered  Fine  enough  for  such  a  calculation  and,  also,  the  outlet  boundary  conditions  are 
imposed  at  a  not  sufFiciently  downstream  location.  A  very  important  point  needs  some  attention: 
the  solutions  of  Ref.  7  use  a  correct  IBL  model  which  contains  the  same  terms  of  the  NS 
equations  considered  by  the  present  IBL-type  approximation,  when  an  appropriate,  orthogonal 
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grid  is  employed.  Therefore,  for  the  present  problem,  the  normal  pressure  gradient  is  indeed 
negligible  in  the  (separated)  boundary  layer,  but  the  present  IBL-type  model  cannot  be  accurate 
insofar  as  it  neglects  the  pressure  gradient  in  the  y  direction,  which  is  considerably  inclined  with 
respect  to  the  direction  normal  to  the  wall,  especially  near  the  x  =  0  and  x  =1  locations.  In 
conclusion,  the  PNS  equations  can  still  be  used  in  the  presence  of  mild  skewness  of  the 
computational  mesh,  but  the  IBL  equations  strictly  require  a  body  oriented  coordinate  mesh 
which  is  orthogonal,  at  least  near  the  body  surface. 

As  far  as  the  efficiency  of  the  calculations  is  concerned,  the  following  points  emerge:  the 
present  approach  is  reasonably  efficient,  but  seems  to  require  from  two  to  four  times  more 
iterations  than  other  current  PNS  or  IBL  solvers  [7-10].  However,  the  present  method  is 
extremely  robust,  insofar  a  no  parameter  was  ever  optimized,  or  changed  when  solving  the  second 
more  difficult  problem,  for  which  the  convergence  rate  was  found  to  be  only  marginally  slower 
than  for  the  first  one  (see  Figures  7  and  9).  Also,  the  present  approach  does  not  resort  ' 
kind  of  inner  iteration,  insofar  as  the  nonlinearities  and  the  ” elliptic”  nature  of  the  flo-  d  are 
accounted  by  a  single  global  iteration  process.  Furthermore,  the  present  method  is  tht  jnly  one 
which  effectively  employs  second-order-accurate  central  differences  for  the  first  derivativ  .  the 
streamwise  direction  and  therefore  is  anticipated  to  be  more  accurate  (for  an  equivalent  mesh)  to 
methods  resorting  to  first-order-accurate  upwind  differences  [8]  or  to  the  FLARE  approximation 
in  the  separated  region  [7].  Finally,  the  convergence  of  the  proposed  approach  can  be  improved 
by  using  backward  sweeping  (at  least  in  the  separated  region)  at  successive  iterations,  horizontal 
line  relaxations  every  once  in  a  while,  or  a  multigrid  approach  [19].  None  of  these  possible 
improvements  was  employed  here,  because  the  aim  of  this  work  was  by  no  means  to  provide  the 
"  most  efficient"  high  Re  weakly  separated  flow  solver,  but  to  demonstrate  that  the  NS  and  PNS 
equations  can  be  solved  with  essentially  the  same  efficiency,  so  that  it  may  be  wise  to  consider  all 
of  the  NS  terms,  in  order  to  avoid  unpleasant  "surprises”.  From  this  respect,  it  is  felt  that  such 
a  goal  has  been  achieved,  together  with  the  secondary  ones  of  providing  separated  flow  solutions 
to  the  classical  BL  equations,  which  are  free  of  the  separation  singularity,  and  of  assessing  the 
influence  of  skewness  in  the  computational  mesh  on  the  performance  of  PNS-type  and  IBL-type 
models. 
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Figure  1.  Channel  flow  geometry  and  mesh  for  Re  «=  10  and  boundary  conditions. 


Figure  2.  Longitudinal  step  size  study  for  Re  =  10  channel  flow. 
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Figure  3.  Transversal  step  size  study  for  Re  •*  10  channel  flow. 
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Figure  6b.  Diffuser  geometry  and  mesh  for  the  case  A  =  -  0.16. 
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Figure  7.  Convergence  history  for  diffuser  flow  (Re  *  6250;  A  * 
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Figure  8.  Wall  vorticity  results  for  diffuser  flow  (Re  *  6250;  A  *  -  0.08). 
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